## vector with package names
x <- c("pbapply", "viridis", "knitr", "kableExtra", "ggplot2", "ape")

aa <- lapply(x, function(y) {
  
  # check if installed, if not then install 
  if (!y %in% installed.packages()[,"Package"]) 
    install.packages(y) 

  # load package
  try(require(y, character.only = T), silent = T)
})
knitr::opts_knit$set(root.dir = normalizePath(".."))

knitr::opts_chunk$set(dpi = 58, fig.width = 18, fig.height = 10, echo = TRUE) 

options(knitr.kable.NA = '')

source('~/Dropbox/Projects/lbh_cultural_evolution/source/custom_treespace.R')

theme_set(theme_classic(base_size = 16))
rb.trees <- readRDS("./output/revbayes_output_in_single_R_object.RDS")

rb.trees <- rb.trees[grep("run_", names(rb.trees), invert = TRUE, value = TRUE)]

trees_diags <- data.frame(model = names(rb.trees))

# get lek name
trees_diags$lek <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 1)

# substitution model
trees_diags$subs <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 2)

# period
trees_diags$period <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 3)

# and which fossils were used
trees_diags$fossils <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 4)

trees_diags$align <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 5)

trees_diags$n_trees <- sapply(rb.trees, function(x) length(x$trees))

# Number of models by iterations and lek
kbl <- knitr::kable(as.matrix(table(trees_diags$n_trees, trees_diags$lek)), caption = "Number of models with a specific number of trees by lek")

kableExtra::kable_styling(kbl)
Number of models with a specific number of trees by lek
BR1 CCE HC1 SUR TR1
3002 12 0 24 0 12
6002 0 24 0 24 0

Tree spaces by leks

All fossils using most recent song types

  • Those printed as NULL did not work (too few tips)
#selected leks
sel_leks <- c("SUR", "CCE", "HC1", "BR1", "TR1")
# sel_leks <- c("SUR", "TR1")

tree_names <- grep("early", names(rb.trees), invert = TRUE, value = TRUE)

pnts_lks <- pblapply(sel_leks, cl = 1, function(i) {
  
  lek_trees <- grep(i, tree_names, value = TRUE)
  # lek_files <- rb.trees[lek_files]
  
  multi_list <- lapply(lek_trees, function(x) read.tree(file.path("./output/most_recent_revbayes_models/", x), keep.multi = TRUE))
  
  names(multi_list) <- sapply(basename(lek_trees), function(x) paste(strsplit(x, "_", fixed = TRUE)[[1]][1:5], collapse = "_"))

  pnts <- try(custom_treespace(multi_list, n.points = 100, method = "RF"), silent = TRUE)

  if (is(pnts, "try-error")) pnts <- NULL
    
  #   
  # gg <- ggplot(data = pnts, aes(x = x, y = y, fill = generation)) + 
  #   geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) + 
  #   scale_colour_gradient(low = "red", high = "yellow") +
  #   geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) + 
  #     theme(panel.background = element_blank(), axis.line = element_line(color = "grey"), 
  #     panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) + 
  #   facet_wrap(~chain,nrow = round(sqrt(length(unique(pnts$chain))))) + 
  #     scale_fill_gradientn(colours = viridis(256)) +
  #    labs(x = "MDS_v1", y = "MDS_v2") else
  #   # ggtitle(label = sprintf("Tree space for %d trees", nrow(pnts) / length(unique(pnts$chain)))) else
  #     gg <- NULL

return(pnts)
  
})

names(pnts_lks) <- sel_leks

saveRDS(pnts_lks, file = "./output/rwty/tree_distance_all_fossils_pairwise_shared.RDS")
readRDS("./output/rwty/rwty_tree_distance_plots_all_fossils.RDS")
## $SUR

## 
## $CCE
## NULL
## 
## $HC1

## 
## $BR1

## 
## $TR1

All fossils using tips shared by 90% of trees

#selected leks
sel_leks <- c("SUR", "CCE", "HC1", "BR1", "TR1")
# sel_leks <- c("SUR", "TR1")

tree_names <- grep("early", names(rb.trees), invert = TRUE, value = TRUE)

ggs <- pblapply(sel_leks, cl = 1, function(i) {
  
  lek_trees <- grep(i, tree_names, value = TRUE)
  # lek_files <- rb.trees[lek_files]
  
  multi_list <- lapply(lek_trees, function(x) read.tree(file.path("./output/most_recent_revbayes_models/", x), keep.multi = TRUE))
  
  names(multi_list) <- sapply(basename(lek_trees), function(x) paste(strsplit(x, "_", fixed = TRUE)[[1]][1:5], collapse = "_"))

  pnts <- try(custom_treespace(multi_list, n.points = 100, method = "shared"), silent = TRUE)

  if (!is(pnts, "try-error"))
  gg <- ggplot(data = pnts, aes(x = x, y = y, fill = generation)) + 
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) + 
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) + 
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"), 
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) + 
    facet_wrap(~chain,nrow = round(sqrt(length(unique(pnts$chain))))) + 
      scale_fill_gradientn(colours = viridis(256)) +
   labs(x = "MDS_v1", y = "MDS_v2") else
    # ggtitle(label = sprintf("Tree space for %d trees", nrow(pnts) / length(unique(pnts$chain)))) else
      gg <- NULL

return(gg)
  
})

names(ggs) <- sel_leks

# saveRDS(ggs,file = "./output/rwty/rwty_tree_distance_plots_all_fossils_shared.RDS")
readRDS("./output/rwty/rwty_tree_distance_plots_all_fossils_shared.RDS")
## $SUR

## 
## $CCE

## 
## $HC1

## 
## $BR1

## 
## $TR1

Early fossils using tips shared by 90% of trees

#selected leks
sel_leks <- c("SUR", "CCE", "HC1", "BR1", "TR1")
# sel_leks <- c("SUR", "TR1")

# tree_names <- grep("early", names(rb.trees), value = TRUE)
tree_names <- names(rb.trees)
pnts_lks <- pblapply(sel_leks, cl = 1, function(i) {
  
  lek_trees <- grep(i, tree_names, value = TRUE)
  # lek_files <- rb.trees[lek_files]
  
  multi_list <- lapply(lek_trees, function(x) read.tree(file.path("./output/most_recent_revbayes_models/", x), keep.multi = TRUE))
  
  names(multi_list) <- sapply(basename(lek_trees), function(x) paste(strsplit(x, "_", fixed = TRUE)[[1]][1:5], collapse = "_"))

  pnts <- try(custom_treespace(multi_list, n.points = 100, method = "RF"), silent = TRUE)

  if (!is(pnts, "try-error"))
  return(pnts) else return(NULL)
#   gg <- if (!is(pnts, "try-error")) ggplot(data = pnts, aes(x = x, y = y, fill = generation)) + 
#     geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) + 
#     scale_colour_gradient(low = "red", high = "yellow") +
#     geom_point(shape = 21, size = 4, colour = "white") + 
#       theme(panel.background = element_blank(), axis.line = element_line(color = "grey"), 
#       panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) + 
#     facet_wrap(~chain, nrow = round(sqrt(length(unique(pnts$chain))))) + 
#       scale_fill_gradientn(colours = viridis(256)) +
#    labs(x = "MDS_v1", y = "MDS_v2") else
#     # ggtitle(label = sprintf("Tree space for %d trees", nrow(pnts) / length(unique(pnts$chain)))) else
#       NULL
# 
# return(gg)
  
})

names(pnts_lks) <- sel_leks

saveRDS(pnts_lks, file = "./output/rwty/tree_distance_fossils_pairwise_shared.RDS")
readRDS("./output/rwty/rwty_tree_distance_plots_early_fossils.RDS")
## $SUR

## 
## $CCE

## 
## $HC1

## 
## $BR1

## 
## $TR1

R session information

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=pt_PT.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_PT.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_PT.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ape_5.4-1         ggplot2_3.3.2     kableExtra_1.3.0  knitr_1.30       
## [5] viridis_0.5.1     viridisLite_0.3.0 pbapply_1.4-3    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5       highr_0.8        pillar_1.4.6     compiler_4.0.3  
##  [5] tools_4.0.3      digest_0.6.27    nlme_3.1-149     lattice_0.20-41 
##  [9] evaluate_0.14    lifecycle_0.2.0  tibble_3.0.4     gtable_0.3.0    
## [13] pkgconfig_2.0.3  rlang_0.4.9      rstudioapi_0.11  yaml_2.2.1      
## [17] parallel_4.0.3   xfun_0.19        gridExtra_2.3    withr_2.3.0     
## [21] stringr_1.4.0    dplyr_1.0.2      httr_1.4.2       xml2_1.3.2      
## [25] generics_0.1.0   vctrs_0.3.6      webshot_0.5.2    grid_4.0.3      
## [29] tidyselect_1.1.0 glue_1.4.2       R6_2.5.0         rmarkdown_2.4   
## [33] farver_2.0.3     purrr_0.3.4      magrittr_2.0.1   scales_1.1.1    
## [37] ellipsis_0.3.1   htmltools_0.5.0  rvest_0.3.6      colorspace_1.4-1
## [41] labeling_0.3     stringi_1.5.3    munsell_0.5.0    crayon_1.3.4